rm(list=ls())
#=====================================================
# Coefficients Plot
#=====================================================

library(arm)
library(foreign)

# set up the working directory in your computer
setwd("")

# =================================================================================
# Figure S1 : period varitions for SEI effects on offspring sex
# =================================================================================

pdf("years_sei.pdf", paper="special", width=7, height=7)
layout(matrix(c(1,2,3,3), 2, 2, byrow=T))


coef.all <- read.csv("year_LPM_sei_main.csv")


coef.t <- subset(coef.all,o_female==1)
coef.t1 <- subset(coef.t, o_dv=="px_daughter")
coef.t2 <- subset(coef.t, o_dv=="fx_daughter")

  d1 <- subset(coef.t1, is.na(o_year)==F)
  d2 <- subset(coef.t2, is.na(o_year)==F)
  
  b0.all <- d1$o_coef;se0.all <- d1$o_se;sig0.all <- ifelse(d1$o_p<0.05,1,ifelse(d1$o_p<0.1,2,0))
  b1.all <- d2$o_coef;se1.all <- d2$o_se;sig1.all <- ifelse(d2$o_p<0.05,1,ifelse(d2$o_p<0.1,2,0))
  
    
  xid0 <- d1$o_year
  xid1 <- d2$o_year

  y0.color <- rep("black",length(b0.all))
  y1.color <- rep("darkgray",length(b1.all))
  y0.lty <- ifelse(sig0.all>0,1,2);  
  y1.lty <- ifelse(sig1.all>0,1,2);  
  y0.pch <- ifelse(sig0.all>0,20,1);
  y1.pch <- ifelse(sig1.all>0,20,1);
plot(xid0, b0.all, xlab = "Survey Year", xlim = c(min(d1$o_year),max(d1$o_year)), 
  main="Among Women : SEI effects on offspring sex", cex.main=1,
     ylim=c(-2,2), ylab="Coef. with 95% CI",cex=1,pch=y0.pch, pty="s",col = y0.color, xaxt="n")
for (i in 1:length(b0.all)) lines(rep(xid0[i],2),
                                  c(b0.all[i]-1.96*se0.all[i],b0.all[i]+1.96*se0.all[i]),
                                  col=y0.color[i], lty=y0.lty[i])
points(xid1+0.3, b1.all, cex=1,pch=y1.pch, pty="s",col = y1.color, xaxt="n")
for (i in 1:length(b1.all)) lines(rep(xid1[i]+0.3,2),
                                  c(b1.all[i]-1.96*se1.all[i],b1.all[i]+1.96*se1.all[i]),
                                  col=y1.color[i], lty=y1.lty[i])

abline(h=0)
axis(1, at=xid0,labels=xid0, col.axis="black", las=2,cex.axis=0.7)


coef.t <- subset(coef.all,o_female==0)

coef.t1 <- subset(coef.t, o_dv=="px_daughter")
coef.t2 <- subset(coef.t, o_dv=="fx_daughter")

  d1 <- subset(coef.t1, is.na(o_year)==F)
  d2 <- subset(coef.t2, is.na(o_year)==F)
  
  b0.all <- d1$o_coef;se0.all <- d1$o_se;sig0.all <- ifelse(d1$o_p<0.05,1,ifelse(d1$o_p<0.1,2,0))
  b1.all <- d2$o_coef;se1.all <- d2$o_se;sig1.all <- ifelse(d2$o_p<0.05,1,ifelse(d2$o_p<0.1,2,0))
  
    
  xid0 <- d1$o_year
  xid1 <- d2$o_year

  y0.color <- rep("black",length(b0.all))
  y1.color <- rep("darkgray",length(b1.all))
  y0.lty <- ifelse(sig0.all>0,1,2);  
  y1.lty <- ifelse(sig1.all>0,1,2);  
  y0.pch <- ifelse(sig0.all>0,20,1);
  y1.pch <- ifelse(sig1.all>0,20,1);
plot(xid0, b0.all, xlab = "Survey Year", xlim = c(min(d1$o_year),max(d1$o_year)), 
  main="Among Men : SEI effects on offspring sex", cex.main=1,
     ylim=c(-2,2), ylab="Coef. with 95% CI",cex=1,pch=y0.pch, pty="s",col = y0.color, xaxt="n")
for (i in 1:length(b0.all)) lines(rep(xid0[i],2),
                                  c(b0.all[i]-1.96*se0.all[i],b0.all[i]+1.96*se0.all[i]),
                                  col=y0.color[i], lty=y0.lty[i])
points(xid1+0.3, b1.all, cex=1,pch=y1.pch, pty="s",col = y1.color, xaxt="n")
for (i in 1:length(b1.all)) lines(rep(xid1[i]+0.3,2),
                                  c(b1.all[i]-1.96*se1.all[i],b1.all[i]+1.96*se1.all[i]),
                                  col=y1.color[i], lty=y1.lty[i])

abline(h=0)
axis(1, at=xid0,labels=xid0, col.axis="black", las=2,cex.axis=0.7)


coef.t <- read.csv("year_LPM_sei.csv")

coef.t1 <- subset(coef.t, o_dv=="fx_daughter")
coef.t2 <- subset(coef.t, o_dv=="px_daughter")

  d1 <- subset(coef.t1, is.na(o_year)==F)
  d2 <- subset(coef.t2, is.na(o_year)==F)
  
  b0.all <- d1$o_coef;se0.all <- d1$o_se;sig0.all <- ifelse(d1$o_p<0.05,1,ifelse(d1$o_p<0.1,2,0))
  b1.all <- d2$o_coef;se1.all <- d2$o_se;sig1.all <- ifelse(d2$o_p<0.05,1,ifelse(d2$o_p<0.1,2,0))
  
    
  xid0 <- d1$o_year
  xid1 <- d2$o_year

  y0.color <- rep("black",length(b0.all))
  y1.color <- rep("darkgray",length(b1.all))

  y0.lty <- ifelse(sig0.all>0,1,2);  
  y1.lty <- ifelse(sig1.all>0,1,2);  
  y0.pch <- ifelse(sig0.all>0,20,1);
  y1.pch <- ifelse(sig1.all>0,20,1);
plot(xid0, b0.all, xlab = "Survey Year", xlim = c(min(d1$o_year),max(d1$o_year)), 
  main="Coefficients on the interaction bewteen SEI scores and parental sex", cex.main=1,
     ylim=c(-2,2), ylab="Coef. with 95% CI",cex=1,pch=y0.pch, pty="s",col = y0.color, xaxt="n")
for (i in 1:length(b0.all)) lines(rep(xid0[i],2),
                                  c(b0.all[i]-1.96*se0.all[i],b0.all[i]+1.96*se0.all[i]),
                                  col=y0.color[i], lty=y0.lty[i])
points(xid1+0.3, b1.all, cex=1,pch=y1.pch, pty="s",col = y1.color, xaxt="n")

for (i in 1:length(b1.all)) lines(rep(xid1[i]+0.3,2),
                                  c(b1.all[i]-1.96*se1.all[i],b1.all[i]+1.96*se1.all[i]),
                                  col=y1.color[i], lty=y1.lty[i])

abline(h=0)
axis(1, at=xid0,labels=xid0, col.axis="black", las=2,cex.axis=0.7)


dev.off()
